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A comprehensive heat and mass transfer computational model is used to analyse the intricate two-way 
coupling arising from the activated chemical reactions involved in the heat treatment of wood. The 2D 
version of a drying code known as TransPore is used to simulate the coupled heat and mass transfer phe¬ 
nomena. This code accounts for the internal pressure in the porous medium. The pyrolysis model describ¬ 
ing the chemical reactions occurring in the main constituents within the cell walls of wood (cellulose, 
hemicelluloses and lignins) is derived using data taken from the literature. Refined computational strat¬ 
egies were required to address the two-way coupling between the heat and mass transfer and chemical 
mechanisms, including the thermal activation of the chemical reactions, together with the treatment of 
heat sources (or sinks) and the production of volatiles. The experimental set-up allows the overall weight 
loss, and the internal temperature and pressure at specific locations within the board to be determined 
during processing. The reported simulations highlight that the model is able to capture two particular 
phenomena observed during the heat treatment of wood: the double pressure peak due to water evap¬ 
oration and volatiles production; and the temperature overshoot during the heat treatment phase. 

© 2009 Elsevier Ltd. All rights reserved. 


1. Introduction 

The motivation for this research arises from the potential usage 
of heat treatment processes to improve some of the characteristics 
of the final wood product, such as its durability and dimensional 
stability [1-3]. In particular, this improved dimensional stability 
and durability can be achieved without the use of external chemi¬ 
cal products, which enables the treated wood to remain an envi¬ 
ronmentally friendly product. The main transport phenomenon 
evident during this process is a thermal flux that induces both 
physical and chemical modifications inside the wood cell structure. 
The rate of pyrolysis is controlled to exploit the thermal degrada¬ 
tion of the wood constituents without the use of oxygen, which 
therefore avoids oxidation. The pyrolysis is carried out in the range 
200-260 °C giving the wood new physical, chemical and mechan¬ 
ical properties. 

The concept of thermal treatment to stabilise the wood 
structure is by no means new. The first pilot studies of the subject 
were reported by Stamm [1] and Buro [4] after which a variety of 
different non-industrial based studies on the thermal treatment of 
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wood were reported [5-9]. Recent efforts on this technology have 
lead to the development of several treatment processes being 
introduced to the European market, especially in Finland, France 
and some other European countries [10-12] and as a direct result 
of this research, some wood improvement processes are now doc¬ 
umented in patent specifications [13-15]. 

Wood is a complex biological material and the thermal decom¬ 
position of the chemical components occurs via a series of chemi¬ 
cal reactions coupled with heat and mass transfer. In most 
publications relating to the heat treatment of wood, reference is of¬ 
ten made to the improved dimensional stability and increased 
resistance to fungi [16]. The heat treatment process can also negate 
changes in the overall wood characteristics [17]. The chemical deg¬ 
radation of wood occurs in the hemicellulose, cellulose and lignin 
components and the changed wood composition results in a lower 
hygroscopicity with a major influence on both dimensional stabil¬ 
ity and durability. It appears that the higher the treatment temper¬ 
ature, the better the durability, stability and biological properties 
of the product become. However, it should be noted that the over¬ 
all mechanical properties of the wood are weakened under such 
conditions. In order to reach a selective degree of depolymerisation 
of the hemicellulose during the treatment process, relatively mild 
conditions can be applied to limit unwanted side reactions, which 
can influence the mechanical properties negatively [18]. Interest is 
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Nomenclature 



A 

(NxN) matrix 

Greek symbols 

b 

vector 

X 

depth scalar (m) 

c 

molar concentration (mol m -3 ) 


Porosity m 3 m -3 

Cp 

specific heat (J kg -1 K 1 ) 

(p 

phase potential 

D 

g 

h 

diffusivity tensor (m 2 s -1 ) 
gravitational constant (m s -2 ) 
intrinsic averaged enthalpy (J kg -1 ) 

o> 

2 

volumetric source term (W m -3 ) 
conserved quantity 

Thermal conductivity W m _1 °C _1 
dynamic viscosity (kg m _1 s -1 ) 
chemical potential (J kg -1 ) 
intrinsic averaged density (kg m -3 ) 
surface tension (N m _1 ) 
volume fraction 

h 

Kap 

Ah w 

Ahi 

A Hr 

averaged enthalpy of bound liquid (J kg -1 ) 
latent heat of evaporation (J kg -1 ) 
differential heat of sorption (J kg -1 ) 
enthalpy of reaction i (J kg -1 ) 
volumetric source term (W m -3 ) 

P 

Pb 

P 

G 

8 

ki 

reaction rate of the raction i (s -1 ) 

CD 

mass fraction 

K 

kinetic transport tensor 
flux expression (kg m -2 s -1 ) 

Superscripts and subscripts 

J 

a 

air 

% 

effective thermal conductivity (W m _1 K 1 ) 

b 

bound 

K 

relative permeability 

c 

C 

capillary 

cellulose 

1 

absolute permeability (m 2 ) 

eff 

effective property 

k m 

mass transfer coefficient (m s -1 ) 

f 

fixed phase 

M 

molar mass (kg mol -1 ) 

fw 

free water 

n 

unit normal 

fsp 

fiber satured point 

P 

pressure (Pa) 

g 

gas phase 

Pvs 

saturated vapour pressure (Pa) 

H 

hemicellulose 

Q 

external heat transfer coefficient (W m -2 °C _1 ) 

i 

general index for vectors and ith for reaction 

R 

gas constant (J mol -1 K 1 ) 

L 

lignin 

S 

volume saturation 

P 

pressure 

T 

temperature (K) 

s 

solid phase 

t 

time (s) 

sat 

saturated state 

V 

velocity (m s _1 ) 

T 

temperature 

W; 

mass fraction of compound i in the gaseous phase 

V 

vapour 

X 

moisture content 

w 

liquid phase 

X 

molar fraction 

(X) 

atmospheric 


often focused on the drying characteristics as well as the increased 
dimensional stability, but never on the mass transfer. 

In the first part of this paper we summarise the selected set of 
pyrolysis reactions (see Rousset et al. [19]) used for the chemical 
reactions occurring in the three constituents of the wood cell walls. 
In the second part of the paper we present the extended formula¬ 
tion of heat and mass transfer in porous media, which accounts for 
the heat and production of volatiles due to the chemical reactions. 
Next, we describe the numerical strategy used to solve this com¬ 
prehensive set of transport equations. 

The last part of the paper describes the refined experimental de¬ 
vice used to collect the detailed information during the heat treat¬ 
ment, including local values of internal pressure and temperature. 
The discussion highlights that the model is able to capture the glo¬ 
bal behaviour of wood during the heat treatment process, includ¬ 
ing two important phenomena: the double pressure peak due to 
water evaporation and the production of volatiles; and the temper¬ 
ature overshoot during the heat treatment phase. 

2. The pyrolysis model 

During the last five decades a number of theoretical and exper¬ 
imental studies have been carried out to gain a better understand¬ 
ing of the complex mechanisms that arise during the thermal 
degradation of biomass. Some of these models have been devel¬ 
oped specifically to analyse and subsequently prevent fire. For 
example, some authors [20,21] have studied the ignition of fire 
in wood, while others [22,23] have investigated flame propagation. 


Fredlund [24] and Yuen [25] have also studied the resistance of 
building materials from fire damage. Pyrolysis models have been 
used in the past to describe, in a global manner, the lesser-known 
mechanisms of wood degradation. Nowadays it is possible to find 
quite complex models that describe wood pyrolysis based on dif¬ 
ferent reactions occurring in parallel, consecutively or competi¬ 
tively [26-28]. 

Numerous mechanisms have been published to provide a ra¬ 
tional explanation of the decomposition of cellulose, hemicellu- 
loses and lignin [29,30]. The interested reader can find a more 
detailed synopsis of these different models in Gronli [31] or Di 
Blasi [32]. One of the primary objectives of this research is to cou¬ 
ple a detailed pyrolysis model that describes the evolution of the 
major chemical constituents in the wood cell wall during thermal 
treatment, with the complex heat and mass transport computa¬ 
tional model known as TransPore [33,34]. While the coupling of 
heat and mass transfer with chemical reaction models is not new 
at high temperatures 400-800 °C [35,36], there is only a limited 
number of models used to investigate this process at low temper¬ 
atures (<280 °C) [37-39]. In most of these models the diffusion of 
water vapour in the gas phase is neglected and the total pressure 
is approximated using only the water vapour pressure. The gases 
from the pyrolysis reaction are generally not taken into account 
in the formulation [40]. 

Generally two different approaches are used in modelling the 
pyrolysis of biomass. The first consists of predicting the total 
behaviour of wood starting from the evolution of its principal com¬ 
ponents [41] according to the following equation: 
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Wood 


k 

—-—► Charcoal 

k 2 


' Tar 


► Gases 


Tar 


—^—►Charcoal 


——►Gases 


Primary reactions Secondary reactions 

Fig. 1. Reactive model for the Lumped Parameter Approach. 


Wood = %cellulose + %lignin + %hemicellulose 


where (%) represents the contribution for each compound in the 
composition of the wood. 

The second method, called the “Lumped Parameter Approach”, 
is more of a global strategy that consists of classifying the products 
resulting from the degradation of biomass by considering that 
there is only one homogeneous element containing condensable 
products (tar), non condensable (gas) and solids (charcoal). The 
model is based on three parallel and competitive reactions, as de¬ 
scribed previously in the literature [42,43]. The presence of sec¬ 
ondary reactions explains the exothermic phenomena [28] (Fig. 1). 

Although this method is the most commonly used, it was not 
retained here. Our goal is to develop a model that would be easily 
applied to other lignocellulosic biomass. This is why we preferred 
to model the thermal decomposition of each component of wood. 
The model adopted here (see Table 1) is based on the in-depth 
analysis presented by Rousset et al. [19]. Observe in Table 1 for cel¬ 
lulose that we have retained the reaction pathway suggested by 
Broido-Shafizadeh [44] and modified by Varhegyi et al. [45], which 
states that the thermal decomposition of cellulose is represented 
by two competitive reactions. Pure microcrystalline cellulose, avi- 
cel PHh-105 was used in these works. To investigate the pyrolitic 
events that occur at lower temperature, an experiment was per¬ 
formed at 5 °C/min up to 370 °C. This simplified model is generally 
well accepted today as the reference point for the majority of work 
relating to the pyrolysis of this compound. The kinetic parameters 
necessary for this reaction have been taken directly from Bradbury 
[46]. Note that the double mechanism reaction is still utilised even 
though in the zone of the temperature of interest the mono-reac¬ 
tion would be sufficient. 

Table 1 also shows an irreversible single order reaction for the 
lignin component, which is generally the accepted model most fre¬ 
quently used to describe its decomposition. By extrapolating the 
losses of mass observed during the experiments [19], we decided 
to use the model proposed by Williams and Besler [47]. In this 
work, the lignin was extracted from pine wood and heated from 
0 °C to 720 °C at 5 °C/min using nitrogen. 


Hemicelluloses are the most thermally unstable compounds and 
are strongly implied in the reactions of degradation and dehydra¬ 
tion. In this work we retained the successive reactions suggested 
by Varhegyi et al. [48] because this model seems more adapted to 
describe the thermal decomposition of hemicelluloses in the range 
200-260 °C. We have compared several models to confirm these re¬ 
sults, however we should not rule out the possibility of adding a 
third reaction, as proposed by Ward and Braslaw [49]. 

From the reaction diagrams given in Table 1 it is possible to 
deduce the following system of differential equations for the 
degradation of the wood components: 

- Hemicellulose component 

^=-mp„, % L=0 - 43 fe i P H - 

%- = -fe ■ Ps, + 0.57 • fc, • p H , ^ = 0.56 • k 2 ■ p Si , 

^ = 0.44 .k 2 ps, (1) 

where p H represents the mass of hemicellulose; kj and k 2 are the 
reaction rates for hemicellulose; and p s . } p G are, respectively, the 
masses of the coal or tar and volatile matters produced during reac¬ 
tion i = 1,2. 

- Cellulose component 

^ = -(k 3 + k 4 )p c , ^ = (2) 

where p c is the mass of cellulose and k 3 and k 4 are the reaction rates 
for cellulose. 


- Lignin component 


dfh 

dt 


“Ml, 


dPc 5 

dt 


Ml, 


( 3 ) 


where p L is the mass of lignin and k 5 is the reaction rate for lignin. 
In the above equations, we model the reaction rates according to 
the Arrenhius Law k, = A, exp (^j where the parameters A* and 
E i} i = 1,.... 5 have been taken from the literature and are summa¬ 
rised in Table 1. These parameters are assumed to be constant in 
the range 180-260 °C. 


3. The heat and mass transfer model 


3.1. Conservation equations of heat and mass transfer 

The macroscopic conservation equations that govern the heat 
and mass transfer phenomena that arise in porous media during 


Table 1 

Parameters used for the wood pyrolysis model. H, C, L designate, respectively, hemicellulose, cellulose, lignin; G, and S, designate the products from the reaction i, respectively, in 
the gaseous phase and in the solid phase. 


Component/authors 

Models 

Kinetics parameters 

Temperature (°C) 

Heating rate (°C min ’) 

Hemicelluloses (Varhegyi) [46] 

H 

1 ^ 

-► 0,43 Gi + 0.57 Si 

2 1 

0,56 G 2 + 0.44 S 2 

F^^kjmor 1 

A 1 = 7, 94 10 16 s -1 

E 2 = 95 kj mol -1 

A 2 = 5,01 10 6 s -1 

250-280 

10 

Cellulose (Varhegyi) [43] 

C 

3 

1- ► S 3 

4 H ► G 4 

E 3 = 147 kj mol -1 

A 3 = 2.51 10 9 s” 1 

E 4 = 238 kj mol” 1 

A 4 = 1.25 10 18 s -1 

259-341 

10 

Lignin (Williams) [45] 


' $5 + G 5 

£ 5 = 124.3 kjmor 1 

A 5 = 2.77 10 7 s” 1 

720 

5 
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drying are now well known [33,50]. This now classical and well ac¬ 
cepted formulation has been extended here to account for the ef¬ 
fects of chemical reactions (component degradation, heat source, 
heat sink and volatiles production): 

- Solid conservation 

p s s ^ + V • (Tvvs) = -£ «m Ci > + <tn Si », (4) 

i 

where (; rh c .) is the volumetric mass rate of gas production, (m s .) is 
the volumetric mass rate of charcoal or tar production, both of 
which are in kg m -3 s -1 and produced during the decomposition 
reaction i. Note that wood deformations are neglected in this work, 
implying that the solid phase velocity v s is assumed zero through¬ 
out the process. 

- Liquid conservation 

Q _ _ _ _ _ _ _ _ — 

g~ t (e w pZ + £gpf, + Pb) + V • (p"v w + pf,Vg + p b \ b ) = V ■ (pfD eJj r' Va> v ). 

(5) 


- Volatile matters and air conservation 

The volatile matters (gases) produced by the chemical reactions 
manifest as source terms for the gaseous phase. The conservation 
equation on the ith gas component G, can be written as 


PqVc, = p g Gl Vg - p|D eff . Vco Gl . 


( 6 ) 


In order to have not more than three dependent variables (tem¬ 
perature, moisture content and gaseous pressure) computed by the 
computational model, we have assumed the mixing gas to be 
homogeneous. This implies that the diffusion of these components 
in air is negligible. The Peclet number is the relevant criterion to 
justify this simplification. In the case of a porous medium, this 
number involves the permeability and the effective mass diffusiv- 
ity. Its value can be estimated by assuming a linear pressure field 
within the sample: 


LV /CAP 1 /CAP 

J 6 “ ~D ~ 1 X ~pL X Df f ~ JiDf f ' 


(7) 


As a result of these assumptions, the air and gases conservation 
law can be written as 


| (%Pf) + V • (Pfv g ) = V ■ (p|D e# V® 0 ) + £ 

° i=1,2,4,5 


M a (ih Ci ) 

M Cl 


( 12 ) 


- Energy conservation 

^ ( e wPw^w + %(Pf hv + Pa fr a ) + Pbhb + P 0 h s — SgPg) 

+ v • (pZh w v w + Q5f h v + p gt h*)v g + h b p,^v~ b ) 

= V • ( plD eff (h v S/(0 v + h a VcOa) + K eJ yVT) - A Hr (13) 


The term A H R within Eq. (13) represents the net energy gener¬ 
ated as a result of the chemical reaction of the constituents of the 
wood cell wall during thermal treatment. Note that the value of the 
enthalpy of reaction i, A h it is assumed constant in the temperature 
range of interest. The precise form of this term is given as follows: 

A Hr = Ah\k\Ptf + A/i 2 /< 2 p Sl + (A /13 /<3 + A/i 4 /< 4 )p c + Ah^ksp^. (14) 

This term is one of the main links between the transfer phe¬ 
nomena and the chemical reactions. It allows the heat produced 
by the chemical reaction to be accounted for in the energy balance. 
The success and validity of the simulations (temperature overshoot 
and induced chemical activation) strongly depends on how well 
the enthalpies of reaction are known. 

In Eqs. (5), (12), (13) the gas and liquid phase velocities are gi¬ 
ven by the Generalised Darcy’s Law: 


K k 

\ e = — ^V(p e ,V(p e = VP e -p e gX7x where ( = w,g (15) 
Pi 

where the quantities § are known as the phase potentials and % is 
the depth scalar. 

In order to close this system, the following constraints and con¬ 
stitutive relations are needed: 


£ g 

Pv 


(j)Sg^ £ W - (j)Syy^ 

MI p* = M I 

M v ’ a M a 


Sg + S w — 1, Pg — P v + P*, 
P w = Pg ~ Pc(5w, T). 


(16) 


Using the values of beech in the radial direction and for a pres¬ 
sure difference of 10 5 Pa, a typical value computed during the gas 
production phase, the Peclet number obtained with this formula 
is about 150. Note that the value would be smaller, but of the same 
order of magnitude in the tangential direction (^75), as the anisot¬ 
ropy ratios are slightly different for permeability and diffusivity. 

It is therefore valid to assume that all gaseous components have 
the same velocity as the air, namely: 

Pc,v c, = p g c y g - (8) 

The equivalent density p g a * produced by the wood decomposi¬ 
tion is defined to keep the additional rule for the partial pressures: 


(9) 


where Pf* is the apparent partial pressure of the mixture of air and 
volatile components given by, 


p?+E 


TfRT 

M Cl 


( 10 ) 


From this relation, one can derive an expression for the equiva¬ 
lent air density: 




( 11 ) 


where M a is the molecular weight of the air and M Cj the molecular 
weight of the ith volatile component. 


The enthalpy-temperature relations are given by the following 
definitions: 

h s = C ps (T-T r ), h; = C pa (T-T R ), 

hv — hpap + Cpv(T — Tr)i h w = C pw (T — Tr ), (17) 

1 f pb 

hb = h w - — / Ah w dp, h b = h w - Ah w . 

Pb Jo 

The mass fractions of vapour and air in the gaseous phase are 
given, respectively, by 

oj v = — and oj q = —. (18) 

Pg Pg 

The bound liquid flux is assumed to be proportional to a gradient in 
the bound liquid moisture content [33] 

p b \ b = (19) 

The hygroscopic behaviour of the product is defined using the 
classical concept of sorption isotherm (Table 3). 

3.2. Boundary and Initial conditions 

There are two types of boundary conditions that need to be dis¬ 
cussed, namely conditions at an external boundary and conditions 
at symmetry planes. The boundary conditions proposed for the 
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external drying surfaces of the sample are assumed to be of the fol¬ 
lowing form [51] 

J w n = k m cM v In 

Je • n = q(T - Too) + h v k m cM v In 

where J w and J e represent the total liquid and energy fluxes through 
the boundary surfaces and x v and T are, respectively, the molar frac¬ 
tion of the gas vapour and the temperature at the exchange surfaces. 
The corresponding variables with the subscript oo indicate quanti¬ 
ties that are characteristic of the external drying conditions. The 
quantities k m and q are, respectively, the mass and heat transfer 
coefficients that are obtainable from classical boundary layer theory. 

The pressure at the external drying surfaces is fixed at the 
atmospheric value: 

Pg = Poo- (21) 

Symmetry planes are introduced to reduce the overall computation 
times and initially the porous medium has some prescribed mois¬ 
ture and temperature distribution, with the pressure being constant 
throughout at the atmospheric value. 

3.3. Computational strategy 

The numerical procedure used to solve the macroscopic scale 
model has been published extensively over the last decade and 



( * } depending on the convergence rate 


Fig. 2. Coupling between the TransPore code and the pyroysis model. 


the interested reader is referred to some of the more recent publi¬ 
cations of the authors for the finer details of the method used 
[33,34]. Briefly, the finite volume method is implemented to dis- 
cretise the conservation laws. Thereafter, an efficient inexact New¬ 
ton method is used to advance the complicated nonlinear system 
that describes the drying process temporally. Flux limiting is used 
for the spatial weighting schemes for all advection/convection 
terms in the equations in order to reduce numerical dispersion of 
the drying fronts. 

The system of differential equations (l)-(3) describing the 
pyrolysis model is integrated in time using an implicit first order 
scheme that is consistent with the method employed for the trans¬ 
port equations. The resulting discrete system analogue can then be 
expressed in the form: 

^(u) = (F h (u) , F Sl (u), Fs 2 (u) , F Cl (u) , F Gi (u) , F c (u), 

xFc 4 (u),F t (u),F G5 (u)) T = 0, (22) 

where u = (p H , p S] , p S2 , p Gi , p Gi , p c , p Gi , p L , p Gs ) T . 

This nonlinear system must be solved, at each time step, in or¬ 
der to advance all of the pyrolysis variables in time. The solution 
procedure is carried out in two distinct stages via what are known 
as outer and inner iteration. During the outer iteration phase, the 
nonlinear system of equations is linearised according to a Newton 
scheme where the estimate of the solution vector at the (n + l)th 
level is computed from the current solution at the nth level as 
u (n+1) = u (n) + <5u (n) . The linearised system J(u (n) )<5u (n) = -^(u^) is 
then solved for the correction vector <5u (n) . Here J e IR 9x9 represents 
the Jacobian matrix and this system is solved using LU decomposi¬ 
tion with scaled partial pivoting to ensure accurate solutions. 

The coupling strategy between the pyrolysis and transport 
models is devised in order to allow the simulations to proceed at 
a reasonable time scale (Fig. 2). At first the temperature of the 
wood is estimated using TransPore for each control volume and this 
temperature then is used to compute the reaction rates k» i = 
1,..., 5 for each control volume. Next, the pyrolysis model is 
advanced in time for each control volume to produce the heat of 
reaction term A H R (14) and to calculate the gas production/con¬ 
sumption rate, which is subsequently used to advance the heat 
and mass transport model to the next time level. Inner iterations 
are used to ensure convergence of the overall process for a given 
time step. The latter has a value automatically updated according 
to the convergence rate of the inner iterations. These variations 
are allowed within the bounded interval [0.01 ... 10 s]. This decou¬ 
pling strategy balances computation time and accuracy and 


1 -x p 
1 Xp J 




Fig. 3. Computational simulation: the mesh and boundary conditions used in the two- dimensional configuration. 
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Fig. 4. Experimental reactor. 



Thermocouple 


Connector T/P 


Capillary filled 
with oil of silicone 


Coated film 
in silicone 


Longitudinal 

direction 


Metallic parts 
strongly tied 


Fig. 5. Sensor for simultaneous measurement of pressure and temperature at the same location. Note the coated device at the endpieces to avoid any leakage due to 
overpressure. 


appears to work well for all of the tests reported in Section 4. The 
typical 2D computational domain is represented in Fig. 3. 

This model allows the drying and heat treatment simulation of 
one board section consisting of 21 x 41 nodes to be completed in 


less than 1 min on a 2.4 GHz Centrino® processor. We also tested 
the effect of the mesh refinement. For the configurations simulated 
in this paper, the time step is mostly governed by the coupling be¬ 
tween transfers and chemical reactions. As the mesh size has little 
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Fig. 6. Positions of the thermocouples and “pressure-temperature” sensors (units in mm). 


Table 2 

Enthalpy for decomposition reactions of hemicelluloses (reactions 1 and 2), Cellulose 
(reactions 3 and 4) and lignin (reaction 5) used in this work for temperatures below 
260 °C. Positive value indicates endothermic reaction, and negative value exothermic 
reaction. 


Reaction, 

Aft, (kj kg 

1 

+42 

2 

Identified 

3 

+150 

4 

+150 

5 

-233 


effect on the time step, we found that the total computational time 
is approximately proportional to the number of nodes. In terms of 
result accuracy, we found that both the temperature profiles and 
the pressure values were almost unchanged by changing the total 
node number by a factor 2 or 0.5 (less than 0.15 °C of temperature 
shift at core, 0.5 °C at surface and same value of pressure peak at 
core, but shifted by about 1 min,). Significant effects appear clearly 
when the total number of nodes is divided by 4 (times 0.5 in both 
directions). In this case, the temperature difference at surface at¬ 
tains 6 °C during the transient thermal periods, namely at the 
end of the drying plateau at 110 °C. 

3.4. Physical properties 

Most of the key parameters within the macroscopic conserva¬ 
tion Eqs. (5), (12), (13) are strongly dependent on the material 
structure. These properties, together with their postulated func¬ 
tional dependencies, have been taken directly from the literature. 
The reader can refer to the comprehensive list of correlations 
supplied in Perre and Turner [33] or Agoua [52] for the capillary 
pressure; relative permeability in the radial, transverse and longi¬ 
tudinal directions; vapour diffusivity; and thermal conductivity. 
The permeability values come from [53]. 


4. Validation of the coupled model 

4.1. Experimental study 

Fig. 4 presents the configuration used to perform the thermal 
treatment experiments for beech. Beech wood ( Fagus sylvatica L.) 
has been chosen for this study because it is widely used for furni¬ 
ture framing, flooring and engineering purposes, and an additional 
thermal treatment on this species is common in wood the industry. 

This configuration consists of an oven heated by electric heaters 
and an oxygen analyser, both of which are controlled by a com¬ 
puter that monitors the experiment. The heat exchange coefficient 
was experimentally determined in turbulent conditions from a 
previous work [54] as q = 28 W m 2 °C. The endpieces of the board 
are coated and the coating is maintained by two metallic parts 
strongly tied to avoid any leakage due to overpressure (Fig. 5). 
Thermocouples were used to measure the temperature at different 
locations within the sample as shown in Fig. 6. 

Pressure measurements are carried out with a KULITE sensor 
(Series XTC-76A-190-35BARSG) (0-35 absolute bars). The dual 
connector “pressure-temperature” is an improved version of the 
connector originally designed by Perre [51,55], refer to Fig. 5. It 
should be noted that the accurate measurement of both pressure 
and temperature at several locations within the sample is a partic¬ 
ularly delicate operation. Our device allowed these measurements 
to be successfully achieved during the experiment. 

4.2. Comparison of experimental and simulation results 

The reaction enthalpy values are essential to correctly predict 
the temperature field, hence the thermal degradations of a macro¬ 
scopic sample. However, the heats of reaction are often reported as 
being either endothermic, exothermic or a combination of both. 
This lack of precision is often said to be due to the experimental 
conditions, which includes the sample size, the nature of the atmo¬ 
sphere, and the presence of impurities. The published reaction 
enthalpies range from -2300 to 418 kj/kg for wood, from -510 
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Fig. 7. Predicted temperature (a) and relative pressure (b) at the core for different 
enthalpy values of hemicelluloses decomposition. 


to 322 kj/kg for cellulose, from -455 to 79 kj/kg for lignin and from 
-700 to 42 kj/kg for hemicelluloses [19]. With these large varia¬ 
tions, ranging from endothermic to exothermic, it is quite difficult 


to supply the code with relevant parameters. Even the most recent 
works are not able to provide new insights into this question [56]. 
However, in the range of temperature concerned by our work, the 
chemical modifications affect primarily hemicelluloses, the most 
thermally unstable compound. Consequently, we allowed the en¬ 
thalpy of the secondary reaction of the hemicelluloses decomposi¬ 
tion to be tuned for this component and fixed the enthalpy reaction 
of cellulose at 150 kj/kg and the enthalpy reaction of lignin at 
-233 kj/kg, which are values generally accepted in the literature 
(summarised in Table 2). 

Two tests were used in this work to compare simulation and 
experiment. In Test 1, a quartersawn beech board with dimensions 
50 x 150 x 250 mm 3 and an initial MC of 15% is dried at 110 °C 
during 1 h, and then heated at 250 °C during 3 h. In this test, the 
core temperature is used to tune the degree of freedom allowed 
in our model (the enthalpy of the secondary reaction of hemicellu¬ 
loses decomposition). 

In Test 2, a quatersawn beech board with half the thickness of 
sample 1 (24 x 75 x 250 mm 3 ) and an initial MC of 5% is directly 
heated at 250 °C during 100 min. This second test is used to vali¬ 
date the model, namely in the evolution of pressure at the core 
of the sample. The endpieces were coated for both tests, which jus¬ 
tifies the use of the 2D version of TransPore in the board section (ra¬ 
dial/tangential plane). 

Fig. 7 depicts the evolution of the temperature and the dimen¬ 
sionless pressure at the core of a board for different values of the 
enthalpy of the secondary reaction of hemicelluloses decomposi¬ 
tion (-100, -500 and -1000 kj kg -1 ). The input values used in 
TransPore are summarised in Table 3. The board section is 
50 x 150 mm 2 and the initial MC 15%. Fig. 7a exhibits the dramatic 
effect of this reaction enthapy on the temperature evolution. The 
rate of temperature rise and the maximum temperature overshoot 
increases with the enthalpy value (a minus sign indicates an exo¬ 
thermic reaction). As expected, the reaction heat does not affect 
the first pressure peak, due to water evaporation. On the contrary, 
the intensity of the second pressure peak, due to the production of 
gaseous by-products, is strongly affected by the reaction heat, 
through the temperature levelling that activates the reaction rate 
(Fig. 7b). 

Confident with the huge effect of our degree of freedom on the 
simulated fields in our range of temperature, the real climatic con¬ 
ditions recorded during Test 1 was used as input data in the model. 
Being aware of the limited accuracy of the other enthalpy values, 


Table 3 

Input values used in the code TransPore. 



Value 


Unit 

Heat transfer coefficient 

28 


W ITT 2 °C“ 1 

Mass transfer coefficient 

28 x 1(T 3 


m s -1 

Initial temperature 

20 


°C 

Initial moisture content 

5 or 15 


% 

Porosity 

0.6 


(-) 

Over dry density 

600 


kg ITT 3 

Cellulose proportion 

50 


% 

Hemicellulose proportion 

25 


% 

Lignin proportion 

25 


% 


Radial 

Tangential 


Liquid permeability 

12xl0 -16 

3xl0 -16 

m 2 

Gaseous permeability 

6xl0 -16 

1.5xl0“ 16 

m 2 

Liquid relative permeability 

Ki=s% 


(-) 

Gaseous relative permeability 

krg = 1 + (2 ■ Sfw — 3) X Sfw 


(-) 

Bound water diffusion 

D b = -9.9 - 9.8X fa - 4300/7 

D b /2 

m 2 s' 1 

Gaseous diffusion 

D e ff = 1< rg ■ D v 10 3 

Deff 12 

m 2 s -1 

Diffusion of vapour in air 

D v = 2.26 x 10“ 5 ■ (T/273) 1 - 81 ■ P a tmlP 


m 2 s -1 

Sorption isotherm 

P V P VS = 1 - exp(-0.764274 - 3.6787 ■ (. X b /X fsp ) 2 ) 


(-) 

Thermal conductivity 

A eff = 0.14 + 0.3 X 


W m 1 °C“ 1 

Differential heat of sorption 

Ah w =0A-h vap -[(X fsp -X b )IX fsp ] 2 


Jkg- 1 
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Fig. 8. Comparison between experimental and simulated temperature curves at the 
core and beneath the exchange surface of the board for Test 1. Simulation done with 
-850 kj kg -1 for the enthalpy of hemicelluloses decomposition (reaction 2, Table 1 ). 
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Fig. 10. Predicted temperature (a) and relative pressure (b) at the core of different 
board sections. The simulations made with the one-dimensional version of 
TransPore are depicted with a dashed line. Simulation done with -850 kj kg -1 for 
the enthalpy of hemicelluloses decomposition (reaction 2, Table 1). 



Fig. 9. Comparison between experimental and simulated temperature curves (a) 
and dimensionless pressure (b) at the core and beneath the exchange surface for 
Test 2. Simulation done with -850 kj kg -1 for the enthalpy of hemicelluloses 
decomposition (reaction 2, Table 1). 

the reaction enthalpy was not fitted by an inverse method. Instead, 
we just tested different values. A good compromise was obtained 
for an enthalpy of -850 kj kg -1 for the secondary reaction of hemi¬ 


celluloses decomposition. In spite of this adjustment procedure, 
this value is not so far from the value proposed by Bonhke 
(-700 kj kg -1 ) [37]. Using this value, the predicted and measured 
temperatures around 5 mm beneath the surface and the core tem¬ 
perature are depicted in Fig. 8. The evolution shape of both temper¬ 
ature curves is quite well captured by the model. 

The same product parameters were used to simulate Test 2. We 
just changed the geometrical characteristics and the measured gas 
temperature used as boundary conditions to perform this simula¬ 
tion. The temperature evolution is quite similar for the experimen¬ 
tal and simulated results (Fig. 9). Most importantly, we note that 
the simulated temperatures and pressures within the board attain 
approximately the same maximum values as those observed dur¬ 
ing the experiment. An important observation from the experi¬ 
mental pressure evolution is the existence of two distinct peaks: 
the first occurring at 30 h due to bound water evaporation; and 
the second at 75 h due to the production of volatiles as by-products 
of the chemical reactions. It has to be noted that the model is able 
to capture these trends, which proves its predictive capabilities. 
Observe however that the first pressure peak is smaller in the 
experiment than in the simulation. One must remember that the 
simulated pressure is quite sensitive to the value of permeability 
used for the simulations [54], and that this value has huge varia¬ 
tions in wood. In addition, the formation of surface checking can 
increase the permeability close to the exchange surface signifi- 
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Fig. 11. Predicted temperature (a) and relative pressure (b) at the core of different 
board thickness. Simulation done with -850 kj kg -1 for the enthalpy of hemicel- 
luloses decomposition (reaction 2, Table 1). 

candy. These checks open during the drying phase, but are closed 
after due to stress reversal [57]. 

In conclusion, our model is able to capture the behaviour of 
wood during heat treatment, including two features very specific 
to this process: a second pressure peak due to the production of 
volatiles and the temperature overshoot (compared to the external 
temperature) due to the heat of reactions. 

As an example of the model potentials, Fig. 10 depicts a compar¬ 
ison of the core temperature and pressures obtained for 50 mm 
thick boards of different widths, ranging from a square section to 
an infinite width (the ID configuration). These simulations showed 
that the effect of the edges disappears completely for the core tem¬ 
perature and the core pressure only for a width 5 or 6 times larger 
than the thickness. It is interesting to notice that two opposite ef¬ 
fects explain why both the temperature overshoot and the pres¬ 
sure peak depend little on the board width. For a narrow section, 
the presence of lateral exchange faces reduces the process time: 
shorter drying phase, shorter heating phase and shorter treatment 
regime. The gas is produced inside the product over a short time 
interval, that would produce a higher peak, but this gas source 
can escape in two dimensions. For large boards, the gas production 
lasts longer, but the 2D effect for Darcy’s regime becomes weaker 
and eventually vanishes for an infinite width. For the infinite 
width, a third explanation is required to explain why the pressure 
peak is smaller in spite of a pure ID convective flow: the temper¬ 


ature rise is slow enough for one part of gas production to occur 
before the temperature plateau. 

Fig. 11 depicts a comparison of the core values of temperature 
and pressuree obtained for 150 mm width boards of different 
thicknesses. It is interesting to notice that the temperature over¬ 
shoot and the pressure peak now depend strongly on the board 
thickness. The lower heat resistance between core and surface of 
thinner boards limit the temperature peak. In parallel the gas can 
escape easier from thinner boards. The previous results therefore 
confirm the key role of the particle size in the control of heat treat¬ 
ment and promote such a predictive model as an efficient tool to 
tune the control of the industrial process. 

5. Conclusion 

The experimental and theoretical investigations on the thermal 
treatment of wood undertaken in this research have led to the 
development of a computational model that couples the thermal 
decomposition of the cellulose, hemicelluloses and lignin wood 
components with the heat and mass transfer phenomena evolving 
in a board up to about 300 °C. We have exhibited the thermal treat¬ 
ment kinetic curves obtained from a sophisticated experimental 
device that enables the internal temperature and pressure at spe¬ 
cific locations within the board, along with its overall weight loss 
to be measured throughout processing. The comparison between 
experimental and simulated results showed a reasonably consis¬ 
tent agreement between theory and practice. In particular, two 
specific features of heat treatment are nicely captured by the code: 
the second pressure peak due to the production of volatiles and the 
temperature overshoot due to the exothermic reactions. 

Nevertheless, this paper also pointed out that important param¬ 
eters involved in the degradation formulation, namely the reaction 
enthalpies, remain poorly-known. There is therefore a deep need 
for further experimental characterisation of these parameters at 
the micro-particle level. Such a task is part of our projects. 

In closing the research, we believe that this new computational 
tool has great potential to facilitate a better understanding of the 
thermal treatment process and may be used in the future to im¬ 
prove and perhaps optimise this promising industrial process, 
either for material, energy or renewable carbon purposes. Future 
studies will investigate the effect that the strong anisotropy in 
the longitudinal direction has on the thermal treatment simula¬ 
tions, which motivates the development of a completely three- 
dimensional coupled model. At the industrial level, it may also 
be important to model the thermal treatment of an entire stack 
[58,59]. 
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